3*7 * 6 / 



N92-11735 \ 


Nested Ocean Models: Work In Progress 

A. Louise Perkins 

Massachusetts Institute of Technology 




& 


v 


Abstract 

This paper details the ongoing work of combining three existing software programs 
into a nested grid oceanography model. The HYPER domain decomposition program 
[8], the SPEM ocean modeling program [5] and a Quasi- Geostrophy model written in 
England are being combined into a general ocean modeling facility. This facility will 
be used to test the viability and the capability of two-way nested grids in the North 
Atlantic. 


1 Introduction 

We are beginning work on a basin wide coarse grid overlaid with finer grids that follow 
major mesoscale and dynamic features in the North Atlantic basin. The grid management 
will be handled by the HYPER domain decomposition program [8], We will consider several 
combinations of solution methods to be used including nesting a primitive equation fine mesh 
solution method within another primitive equation coarser grid solution, and a primitive 
equation fine mesh solution within a coarser quasi-geostrophic model solution. 

It is well known that to refine the entire coarse mesh in space for ocean circulation 
modeling would be inefficient; it would require large amounts of memory and waste processor 
time in quasi-geostrophic regions. In short, refining the entire coarse mesh is overkill. For 
explicit time-evolution solution methods of primitive equations the advancement must also 
be severely refined in time to account for the gravity wave stability constraint. This results 
in an excessive number of time steps. Alternatively, an implicit solution on a fully refined 
mesh results in a very large matrix problem. 

We are attacking two areas of fundamental ocean modeling directly. The efficiency of 
the boundary conditions between quasi-geostrophic and primitive equation models should be 
advanced based on the insight acquired from our hierarchical approach to the nesting exper- 
iments. The second fundamental area is ocean modeling in general. Nested basin/regional 
grids arc a new concept for ocean applications, and in this respect oceanographic modeling 
lags behind atmospheric and aerodynamic modeling. But the success of domain decomposi- 
tion in these more advanced modeling areas provides encouragement that our research efforts 
are timely, central and well directed towards new, successful applications. 


2 The Navier- Stokes Equations on a Rotating Sphere 

To examine the Navier-Stokes equations on a rotating sphere in a rotating reference frame 
R, let I denote an inertial reference frame, and let r be a radius of that sphere. Then 

(r*)|/ = { r t)\n + ft x r 
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where r t ~ u is velocity and the second term on the right-hand side of the equation is the 
motion a non-rotating observer would see because of the rotation of the sphere. Then 

~ 4~ 2(1 x ur -f" x (f 1 x r) + x v . 

On the right-hand side of this equation the second term is the Coriolis acceleration, the third 
term is the centripetal acceleration, and the fourth term is the acceleration resulting from 
any changes in the rotation speed. 

For geophysical applications here on Earth this last acceleration term is discarded except 
for very long time scales, and centripetal acceleration can be expressed as a potential, 

fix(fixr) = — V0 C , 

that then can be added to the gravitational force potential to net a new geophysical force 
potential. 

The total or material derivative of a scalar quantity is the same in both reference frames. 
Hence the form of the conservation of mass and the thermodynamic equations remains the 
same. 

To estimate the frictional forces F> we could assume a Newtonian fluid with a symmetric 
Navier-Stokes internal pressure tensor. But this molecular dissipative strength would have 
an unknown relationship to the dissipative strength of a given mesoscale ocean phenom- 
ena. In general, a qualitative description of the transfer of energy and momentum between 
scales of interest, and not these smaller molecular scales, are parameterized based on known 
qualitative ocean behavior. 


3 Governing Equations 

« 

Coupling hydrodynamics and thermodynamics, consider an adiabatic, inviscid fluid. It can 
be described by conservation of momentum, a continuity equation, and an energy equation 
coupled with an equation of state. That is, 

d 

—u + aVP + F = 0 
dt 

d „ 

—a — ttV • u — 0 

dt 

4-P + P 7V ■« = ()' 
dt 

is the specific volume, P is the pressure, and 7 is the ratio 
Coriolis and gravity forces, and + u • V is the total 


where u £ R 3 is the velocity, a 
of specific heats. Here F is the 
time derivative. 
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3.1 Primitive Equations 

The hydrostatic approximation neglects vertical acceleration in the vertical equation of mo- 
tion; 

dP =z g 

dz a 

The Boussinesq approximation replaces density with a zeroth order mean density ev- 
erywhere except where multiplied by gravity. The combined hydrostatic and Boussinesq 
approximations are used to formulate a set of reduced equations known as the primitive 
equations. 

These equations are as follows: 

juji + aV H P | F Jf = 0 
dP 

“ 7)7 + 5 = 0 

d 

—a — ov • u = 0 

dt 

d 

— P + P7V • u = 0 
dt 

where u H = (tii.tij)*, u = (u H ,w)\V H = 


3.2 Geostrophy 

Geostrophy, or geopotential flow, retains only the balance between the Coriolis force and the 
potential field. Geostrophy: 

fv = gvx 
fu = gVy 

is the first approximation in an asymptotic expansion of the primitive equations. Here rj is 
the variation of the sea surface height, a measure of pressure. 

3.3 Quasi-Geostrophy 

Large-scale ocean movement is typically quasi-geostrophic. Asymptotically, quasi-geostrophic 
motion has time scales not smaller than the advective time scale. It is 'geostrophic to low- 
est order, yet retains dynamics. Velocity fields can change, but they do so in continuous 
geostrophic balance with pressure. Hence there are temporal derivatives retained. The 
quasi-geostrophic equations are 

u t - fv= -gr) x 
v t + fu = -gg y 

r]t + ( uH) x -|- ( vH)y - 0 . 

Here H + tj is the mean height, //, plus the variation of thickness rj. 
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4 Related Research 


Spall and Holland in unpublished work nested a primitive equation model within ja quasi- 
geostrophic model. They found that the quasi-geostrophic boundary conditions seriously 
dampened the primitive equation physics close to the boundary, reducing them to essentially 
quasi-geostrophic physics in the boundary region. Reducing or eliminating this boundary 
layer is the principle focus of our current efforts. 

Thompson and Schmitz [7] varied the damping time scale on outflow boundary conditions 
for a model of the gulf stream. They found that the outflow dynamics and hence the location 
of the Gulf stream are significantly impacted by the outflow boundary conditions. With 
such strong impact, the possibility of numerical artifacts in regional models due to boundary 
conditions seems large. The lack of existence of well-posed boundary conditions for primative 
equations complicates this problem because no comparisons with true boundary conditions 
can be made even as approximations [7]. See the article by Dr. Blake in this proceedings for 
more detailed information. 

This leads to the importance of getting the boundary conditions physically correct. It is 
known that the subcharacteristics of the Euler equations, upon which the primitive equations 
are based, can have combined inflow and outflow characteristics at both advective inflow and 
outflow boundaries, dependent on the sound speed. Hence the dynamics of the refined regions 
typically has an affect on the surrounding flow fields. Using one-way boundary conditions 
for inflow and one-way for outflow is not, in general, sufficient. There must be a stronger 
interaction between the dynamics of the coarse and the refined meshes. 

To strengthen the interaction Spall and Holland [9] added a direct, but averaged, insertion 
of the streamfunction field (generated using the refined primitive equation depth- averaged 
horizontal velocity values on the refined mesh) onto the coarser quasi-geostrophic solution. 
They allow the refined mesh to dictate the regional external flow component. 

Our hypothesis is that the external flow is insufficient. The strong dynamics in the gulf 
stream are forced by internal instabilities. We are currently testing both baroclinic and 
barotropic nudging for all prognostic variables to establish a stronger relationship between 
the regional dynamics and the coarse mesh solution behavior. This will be done while 
maintaining Spall and Hollands quasi-geostrophic to primitive equation boundary conditions 
[9], which are barotropic, to see if we can improve upon their results. If we are successful, 
we will experiment with a quasi-geostrophic coarse mesh solution and try to find more 
comprehensive two-way communication techniques that reduce or eliminate the boundary 
layer formation found in the previous work. 

We are currently using the barotropic modon defined in [9], and a baroclinic vortex 
problem. 

We are working with flat topography at first, so that we can combine either fixed or 
sigma (stretched) vertical coordinate models. 
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5 Initial Boundary Value Problems 


For modeling purposes these equations must, of course, be viewed as initial boundary- value 
problems (IBVP). Computational models use IB VP for liindcasts, nowcasts, and forecasts 
as well as for physical studies that examine phenomena of interest such as energy cascades, 
eddy shedding, and coastal upwelling. There are two prevalent boundary conditions used in 
ocean models: rigid walls and open boundary conditions. The trick, of course, is to find open 
boundary conditions that are well posed and yet not overconstrained. Because one cannot 
simultaneously diagonalize the coefficient matrices for the multidimensional advective terms 
of the Navier-Stokes equations, this is often a time-consuming guessing game. 

To define open boundary conditions in general, define a system of equations 

Lu = F 


with given initial conditions 

u(x : t = 0) = u,o(^) 
and characteristic boundary conditions 

u l — Su 2 + g 


where u 

[4]- 


(u l ,u 2 ) 1 are the prognostic variables, and S is a generalized reflection operator 


But for reduced equations there can also be modeling constraints such as hydrostatic 
balance or incompressibility. The modeling constraints assumed in order to reduce the 
equations (that is, the asymptotic balances chosen) must be enforced on the initial conditions 
and the boundary conditions to avoid introducing inappropriate length and time scales. 

Many obvious well-posed boundary conditions are overspecified, which leads to the for- 
mation of boundary layers within which the solution adjusts to the additional information. 


5.1 Open Boundary Conditions 

The major issue to address is boundary conditions, Oliger and Sundstrom in [4] detail some 
boundary treatment for geophysical problems, and show that point- wise local boundary con- 
ditions for the primitive equations are not well- posed: the regional open boundary problem 
is open-ended. It is not known, however, whether non-local boundary conditions, such as 
those generated with a domain decomposition method where boundary conditions that are 
derived from a larger domain are or are not well-posed. We may have fewer problems with 
open boundary conditions at a two-sided boundary. The quasi-geostrophic boundary layer 
within the nested primitive equation model in the unpublished work of Spall and Holland 
indicate that open boundary conditions on a nested grid is a major problem that must be 
addressed before two-way nesting will be successful. 

Our approach has been to work from the primitive equation homogeneous model back- 
wards to the quasi-geostrophic model, assessing the differences in the information transmit- 
ted across boundaries, to derive better open boundary conditions for the simplier quasi- 
geostrophic model. 
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For the heterogeneous boundary conditions between the primitive equation and quasi- 
geostrophic regions the quasi-geostrophic boundaries need to evolve as if there was primitive 
equation physics in the region surrounding the refined domain. To insure this we are moni- 
toring §f on the boundary, where quasi-geostrophic physics, which is less vertically diverse, is 
statically stable, and comparing the evolving quasi-geostrophic boundary against a fine mesh 
primitive equation global solution. This is one of our measures of error that is physically 
based. 


0 HYPER 

The HYPER program, described in Perkins [5], looks at domain decomposition as a tool to 
combine grids for computational efficiency and for model flexibility. It currently can locate 
where refined grids should be placed based on asymptotic and-or physical criteria and it 
initializes the refined grids using local uniform mesh refinement. 

Our current work is a static domain decomposition; we are running experiments on the 
influence of the internal boundaries on the flow pattern, and are interested in the flow through 
that boundary. These results are currently being prepared and will appear in a later report. 

The goal is to follow different asymptotic regimes within the ocean basin that are iden- 
tifiable as distinct physical regions of the ocean. The problem is that, once you’re inside 
a reduced physics region, such as a quasi-geostrophic region where there is no ageostrophic 
flow, there may be no way for the model to evolve the complete physics you hope to recover 
by using the refined meshes. For example, in the gulf stream region meanders pinch off to 
form eddies. Many aspects of the physics contribute to this pinching off. In such cases the 
reduced quasi-geostrophic model will not reproduce the ageostrophic behavior in the initial 
conditions of the refined mesh, and hence will miss some of the time dependent interac- 
tions that contribute to the dynamically significant event of ring shedding. There will be no 
ageostrophic signals in the initial conditions to interpolate onto the refined mesh. In such a 
situation, the data used to help initialize the model may make up for some of the missing 
physics, but the time scales may be off. 

7 Domain Decomposition 

Domain decomposition allows the mesh to evolve with the solution. It has been applied to 
elliptic and hyperbolic equations for several years. The interested reader can refer to Chan, 
et al., [2] for an extensive bibliography on elliptic and hyperbolic domain decomposition 
methods. 

Formally we describe the domain decomposition of a discrete coarse mesh, ft 1 , on the 
computational mesh fi, for a fixed time t = nAt, by letting ( p(t ) - 1) be the number of 
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refined subdomains used at time t and {tt k } p ^\ be those subdomains: 

p(0 

U n*(i) eft. 

k~l 

After the domain is decomposed, local uniform spatial mesh refinement, as developed by 
Berger [1], is applied to the new subdomains. The time step may also be refined on these 
regions. 

The sequencing for one coarse time step of magnitude A c t from time t n to time t n + A c t, 
where n indexes the discrete time steps on the coarse mesh, is presented next. Let the 
temporal refinement ratio from the coarse mesh to the refined meshes be r, and notate this 
rAft = A c tj so that a subscript “c” informs us that we are on the coarse mesh, and a 
subscript “/” informs us that we are discussing one of the refined meshes. The domain 
decomposition algorithm follows: 

Domain Decomposition Algorithm 


Advance coarse mesh 

Hark points with signif icant mesoscale and ageostrophic potential 
Cluster these points into refined meshes 
DO r times 

Solve equations on refined meshes 
ENDDO 

Nudge refined values onto coarse mesh 


When all of the refined meshes have been advanced r refined time steps to the next coarse 
time step, their values at time are passed to the coarse mesh where a nudging technique 
modifies the coarse advanced values and produces an aggregate solution on the coarse mesh. 
Let F A , F 1 , and represent the discrete operators for the aggregate solution on the 

original discretized mesh ft 1 , the solution on the coarse mesh ft 1 , and the separate solutions 
on each of the refined subdomains {ft fe }j^, respectively. Then the aggregate solution on the 
coarse discretized mesh is given by * 

F' , (n>) = c[{F‘(n‘)}£ , ,,F 1 (si 1 )] • (7.1) 

where the operator C is a nudging technique that may vary between experiments. 

We use domain decomposition as a tool to combine the explicit coarse mesh solution 
method with the refined mesh solution methods to satisfy our varying numerical requirements 
in a computationally efficient way. We use a two-level refinement scheme consisting of one 
coarse mesh and a set of overlaid refined meshes, where the coarse mesh adequately represents 
quasi-geostrophic behavior, while the refined meshes adequately resolves the more physically 
complete primitive equations. 
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Refined to coarse mesh communication (feedback) can take the form of value averaging, 
as in Berger [1] and as in Spall and Holland [6]. 

We are using a nudging data assimilation technique for the initial experiment and we do 
not include explicit conservation enforcement. 

Our domain decomposition work focuses on two-way interactive nested grid communi- 
cations and the development of good internal boundary conditions. We are particularly 
interested in examining heterogeneous open boundary conditions between different dsymp- 
totic regions. Because the major geophysical equations are not known to be well-posed as 
an initial value boundary problem, in general, this issue becomes important. * 


8 Initial Conditions 

Our long range plans are to build a basin wide grid, and overlay it with refined grids about 
regions of ageostrophic dynamic regions. The refined grids will then follow mesoscale or 
planetary scale dynamical features. 

But our current work is much less ambitious. We have constructed a box model and are 
using Spall and Hollands [12] barotropic inodon and baroclinic vortex problems to examine 
the viability and desirability of different communication schemes between the coarse and 
refined meshes. 

A barotropic modon is a coherent, concentric streamfunction. The barotropic flow is the 
primary mode of a quasi-geostrophic equation formulated as a Sturm-Liouville problem (all 
other modes of the Sturm-Liouville problem are referred to as baroclinic). It has an analytic 
solution, it is quasi-geostrophic, and uses an infinite beta plane approximation. The result 
is a coherent depth independent (barotropic) structure that moves at constant speed. 

The baroclinic vortex has no analytic solution, and is defined using a Gaussian pressure 
distribution with maximum geostrophic velocity of 100 cm! sec. The initial velocity fields are 
calculated to be in geostrophic balance with the prescribed Gaussian pressure field. 

Our experiment follows a hierarchical approach. Beginning with a homogeneous domain 
decomposition we use a full, coarse primitive equation model and keep track of the flow 
across the “future” internal boundaries. Then we introduce the nested grid into the same 
problem and analyze any changes in the boundary information flow. This is used as our 
error due to boundary conditions only. This error is measured both in root mean squared 
error and in phase error. Small shifts of mesoscale features are not always bad compared to 
changes in dynamics within those mesoscale features. 

Once we complete our homogeneous studies we will move to a heterogeneous domain 
decomposition with a quasi-geostrophic coarse grid overlaid with a primitive equation refined 
grid. Again we will compare the flow across the internal boundaries. Then we will add a 
feedback loop that uses the nudging data assimilation technique from the refined mesh to 
the coarse mesh and nudge to the true boundary information. This way we can measure the 
improvement due to the nudging feedback loop. 

The quasi-geostrophic coarser model will be forced by the nudging from the refined prim- 
itive equation model. The unforced quasi-geostrophic model is statically stability, but the 
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primitive equation is not. The forced behavior perturbs the reduced equation dynamics, so 
that a time series of its behavior in the refined region will not be statically stability due 
to the feedback interaction. But with known density changes at the boundaries from our 
true solution, we can measure how well the reduced dynamics are being influenced by the 
regional models with their more complete physical models. Another metric is the apparent 
ageostrophic time series behavior in the quasi-gcostrophic coarser model. After nudging, 
the ageostrophic adjustment to the quasi-gcostrophic coarse mesh is calculated, and a time 
series of this difference is the ageostrophic forcing of the quasi-geostrophic model. Where 
this difference is small, there is no need to maintain a refined mesh, so this metric can be 
used to eliminate refined meshes that are no longer needed, but it can not help us locate 
where refined meshes should be placed. 


8.1 Semi-Spectral Primitive Equation Model (SPEM) 

The primitive equation SPEM model of Haidvogel et. al. [3] has prognostic variables for 
horizontal velocity, u and n, and temperature t y and salinity s . It uses the hydrostatic and 
Boussinesq approximation. The resulting equations are advanced on an scattered Arakawa 
“C” grid in the horizontal, while the vertical is spectral, with Chebyshev modes. It has a 
rigid lid approximation at the surface (no variations or “waves” in sea surface height). 

9 Current Summary 

The computational demands of fully three-dimensional global ocean modeling seem to re- 
quire a nested heterogeneous adaptive grid solution. However, the implementation difficulties 
are robust. The need for physically realistic open boundary conditions is already well doc- 
umented, mostly a result of a “grand challenge” issued several years ago. Our experience 
indicates that an equally pressing need is to provide modeling-consistent asymptotically 
nested initial conditions for each new nested grid. 

The scientific aspects of the work are focused on the boundary condition formulation and 
on the two-way grid communication mechanisms under development. 


References 

[1] M.J. Berger, Adaptive Mesh Refinement For Hyperbolic Partial Differential Equations , 
Ph.D. Dissertation, Stanford University, 1982. 

[2] Chan, T.F., R. Glowinski, J. Periaux, and O.B. Widlund, eds, Domain Decomposition 
Methods , SIAM, 1989. 

[3] Haidvogel, Wilkin, and Young, A semi-spectral primitive equation ocean circulation 
model using vertical sigma and orthogonal curvi-linear horizontal coordinates , JCP, 1989, 
in press. 


G8 


[4] J. Oliger and A. Sundstrom, Theoretical and practical aspects of some initial boundary 
value problems in fluid dynamics , SIAM J. Appl. Math. 35, no. 3 (1978) 419-446. 

[5] Perkins, A. L., Parallel Heterogeneous Mesh Refinement For Multidimensional 
Convection- Diffusion Equations Usinq An Euler-Laqranqe Method, -Ph.D. Dissertation, 
UCRL-53950, 1989. 

[6] Spall, M.A., and W.R. Holland, A Nested Primitive Equation Model For Oceanic Ap- 
plications, to appear in J. Phys. Ocean, 

[7] Thompson, J.D., and W.J. Schmitz, A Limited-Area Model of the Gulf Stream: Design, 
Initial Experiments, and Model-Data Intercomparison, J. Phys. Ocean., Vol. 19, No. 6, 
June 1989. 




69 


